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1 Introduction 

Multisensor electro- and magnetoencephalographic 
(EEG and MEG) as well as electro- and magneto- 
cardiographic (ECG and MCG) recordings have been 
proved useful in noninvasively exctracting informa¬ 
tion on bioelectric excitation [2, 6]. Since the body af¬ 
fects, i.e. filters, the electromagnetic signals measured 
outside of the body, the anatomy of the patient need 
to be taken into account, when excitation sites are lo¬ 
calized by solving the inverse problem. Boundary el¬ 
ement method (BEM) is a widely adopted technique, 
where the body is modeled as an inhomogenous vol¬ 
ume conductor presented by surfaces, such as scalp, 
skull and brain in encephalic applications, or thorax, 
lungs and heart surfaces in cardiac applications. The 
anatomy of the patient has to be individually modeled, 
if accurate high-quality inverse solutions are desired. 
Modern medical imaging modalities, such as mag¬ 
netic resonance (MR) imaging and computerized to¬ 
mography (CT), provide detailed three-dimensional 
(3-D) anatomic information. MR imaging has estab¬ 
lished its position as a standard method to retrieve 
anatomic information for bioelectromagnetic model¬ 
ing. Nevertheless, various solutions to the model con¬ 
struction, in the context of bioelectromagnetic prob¬ 
lems, are not widely reported. A few procedures, 
more or less automated and flexible, have been, how¬ 
ever, described [1,9, 10]. 

This work presents a methodology to construct in¬ 
dividualized boundary element (BE) models for bio¬ 
electromagnetic inverse problems. The objective of 
the work was to develop methods and software, which 
construct boundary element models as automatically 
and fast as possible. This work focus on cardiac prob¬ 
lems, but the modeling of the head is also addressed. 

2 Methods 

The bioelectromagnetic inverse problem pose some 
requirements for the construction of BE models. The 


computation of the inverse problem may become 
a time-consuming process for very dense triangula¬ 
tions. In practice, a few hundred triangles per surface 
is typically used in ECG or MCG calculations yield¬ 
ing a computation time of about 30-60 seconds with a 
standard UNIX workstation. In brain applications, the 
surfaces up to a few thousands of triangles are used. 
In theory, the higher the number of triangles the better 
localization results can be achieved. For numerically 
stable solutions, it is advantageous to compose the tri¬ 
angles as equilateral as possible. Moreover, triangu¬ 
lations where the size of triangles is smallest near the 
heart, are often preferred because the electromagnetic 
gradients show there largest spatial variations. How¬ 
ever, many research groups seem to use contour 1 - 
based triangulation methods, in which the size of the 
triangles is difficult to control. Sometimes, it may 
also be desirable to set triangle points on electrode lo¬ 
cations of a multichannel ECG recording. Moreover, 
the triangulations used for bioelectromagnetic inverse 
problems have to be isosurfaces, i.e. each edge of a 
triangle is shared between two neighboring triangles. 
The methodology to construct individualized bound¬ 
ary element models has been divided into three main 
phases: 1) the segmentation of the MR volume, 2) 
the triangulation of the segmented volume, and 3) the 
registration of the model with the sensors of the bio¬ 
electromagnetic measurement system. 

2.1 Segmentation 

Several methods were tested for segmenting thorax, 
lungs, epicardium and ventricles from MR volumes. 
Finally, a novel method based on deformable mod¬ 
els was developed, described in detail in [4]. In 
the method, a realistically shaped 3-D boundary el¬ 
ement template, consisting of triangulated surfaces, 
is matched elastically to edges extracted from an MR 
volume. The deformation of the model is carried out 

l A contour is a non-overlapping connected point set on a 
plane. 



by the free-form deformation (FFD) approach; the 
model inside a deformation grid is deformed by dis¬ 
placing the points of the grid. The deformation is ac¬ 
complished by minimizing an energy function, 

E — ^data T" 1/Emodel- (1) 

The first term Edata describes the match of the model 
with the salient edges of the images, determined ei¬ 
ther by thresholding or using a standard edge de¬ 
tection algorithm. The second component E mo d e i 
presents how much the model has been deformed 
from its original shape during the process; in practice 
the change in the orientation of the surface normals 
are detected. The user can balance the energy terms 
by the parameter 7 . All grid points are moved sequen¬ 
tially toward the minimum using a conjugate gradient 
method. All surfaces are matched simultaneously. 

In addition, the multiresolution and the global-to- 
local approaches are applied. The multiresolution ap¬ 
proach means that the model is first matched to low- 
resolution images, extracted by low-pass filtering and 
subsampling the original MR volume. As the energy 
minimum is found, the matching is continued at a 
higher resolution level. The matching is more robust 
and faster, when the multiresolution approach is uti¬ 
lized. In the global-to-local approach, more degrees 
of freedom are added gradually to the deformation 
process. In practice, the deformation is started us¬ 
ing a 3 x 3 x 3 grid, and the number of grid points 
is increased as the energy minimum is found. The 
global-to-local approach aims to preserve the original 
shape of the model, representing a typical shape of 
the object of interest, in the beginning of the process 
when the model is usually far from the correct result. 

2.2 Triangulation 

The triangulation can be carried out in two ways: 

A) Since a triangulated boundary element template is 
matched to MR data in the segmentation, the final re¬ 
sult is directly a triangulated surface. The problem 
is that the segmentation algorithm considers only the 
geometric match not the requirements posed for the 
triangulation by the bioelectromagnetic inverse prob¬ 
lem. However, the problem can be partly overcome 
as the deformation defined by the segmentation algo¬ 
rithm is also applied to a separate triangulation satis¬ 
fying the requirements. Moreover, a separate condi¬ 
tion can be added into the energy component Emodei 
(Eq. 1) with the aim of preserving the original dis¬ 
tances between the nodes of the model. 


B) If accurate positioning of the nodes is needed in 
the triangulation, a separate triangulation algorithm 
is applied. The method used in this work is de¬ 
scribed in detail in [5]. The triangulation consists 
of four phases: 1) A binary volume is created from 
the segmentation result for each object of interest. 2 ) 
A marching cubes triangulation consisting of a very 
high number of triangles is constructed for the bina¬ 
rized volume. 3) The nodes desired in final triangu¬ 
lation are chosen from the surface of the marching 
cubes triangulation. 4) A Voronoi area 2 is computed 
for each node, and the neighboring Voronoi areas are 
connecting forming a Delaunay triangulation. The 
Voronoi areas are computed using geodesic distances 
on the marching cubes surface. 

2.3 Registration 

The triangulated surfaces are next transformed, 
i.e. registered, to the position corresponding to the 
pose and the position of the patient under the bioelec¬ 
tromagnetic sensors during the measurements. Since 
the bioelectromagnetic recordings do not provide any 
anatomic information on the positioning of the pa¬ 
tient, a registration method based on external markers 
was adopted. 

Nine external markers are used; five locations are 
chosen in head-feet direction centered to the sternum 
and other four in left-right direction, forming a cross¬ 
shaped pattern. The separation between neighboring 
markers is about 5 cm. In practice, we attach a cross¬ 
shaped object consisting of two strips of rubber on 
the skin. The positions of the markers have been in¬ 
dicated on the strips. The locations of the markers are 
determined relative to the bioelectromagnetic sensors 
using a 3-D digitization system (3SPACE, Polhemus 
Inc., Colchester, VT, USA). 

Plastic cross-shaped markers filled with 
lmmol/lMnCl 2 fluid are attached on the skin 
of the patient before MR imaging. Thus far, the 
markers are manually located from MR images using 
a special software. 

The nine markers located both in the bioelectromag¬ 
netic and MRI coordinate systems are registered us¬ 
ing a non-iterative least-squares method [ 8 ]. The 
transformations discussed here are rigid, i.e. only ro¬ 
tations and translations are concerned. 

2 The Voronoi area of the node is the region closest to this 
node. 



3 Results 

3.1 Segmentation 

So far, the method has been utilized to segment about 
40 patient MR volumes of the thorax. Several MR 
volumes have also been segmented using our head 
and heart models. The computation time is about 1—2 
minutes for a 128 x 128 x 128 volumes using a Sun 
UltralO workstation. 

The matching accuracy of the automated algorithm 
depends on several factors, such as the weight of the 
preservation of the original shape in the minimization 
of the energy, and the final number of grid points, 
which is normally about 15xl5xl5in our stud¬ 
ies. The algorithm produces a good overall match, but 
a few interactive corrections are usually needed. We 
have developed an interactive software based on the 
3-D deformation of the surface at the region of inter¬ 
est. The user-interface is visualized in Fig. 1. In addi¬ 
tion to the interactive software, we have developed a 
modified 3-D version of the original snake-based ap¬ 
proach [3] to fine-tune the result, and to match small 
geometric details, if needed. 



Figure 1: The user-interface dedicated to correct the 
segmentation result. The model superimposed on the 
MR volume is visualized by solid white lines. The 
white 3-D box superimposed on the cross-sections 
shows the region of interest. The deformation of the 
model is accomplished by displacing the center of the 
box; the deformation is equal to the displacement in 
the center of the box and zero on the border of the box. 
In addition to elastic deformations, the model can be 
also traslated, rotated and scaled. 


3.2 Triangulation 

Fig. 2.a shows the triangulation of the thorax pro¬ 
duced directly from the segmentation without an ex¬ 
plicit triangulation method. The triangulation fulfills 
the requirements posed by the bioelectromagnetic in¬ 
verse problem reasonably well, such as the triangles 
are fairly equilateral and the size of the triangles is 
smallest close to the heart where the electromagnetic 
gradients are highest. 

In Fig. 2.b, the triangulation of the scalp and skull 
are visualized. The triangulation is created using the 
Voronoi-Delaunay- based algorithm. The computa¬ 
tion time to build a triangulation from 128x128x128 
binary volume is about 3 seconds. Therefore, the ap¬ 
pliance of the explicit triangulation method does not 
increase the total computation time of the model con¬ 
struction process considerably in practice but gives, 
however, a more optimal triangulation. 

a) b) 


Figure 2: a) The triangulation of the thorax directly 
from the segmentation, b) The boundary element 
model of the scalp and the skull using the Voronoi- 
Delaunay triangulation. 


3.3 Registration 

The user-interface for detecting the cross-shaped MR 
markers from MR images and for registering the 
model with bioelectromagnetic sensors is shown in 
Fig. 3. The root mean square (RMS) error of nine 
registered markers has been about 6 mm in our pa¬ 
tient studies; the lowest value has been about 3 mm 
and the highest more than ten millimeters. 

3.4 Methodology in practice 

The methodology has been successfully applied in 
more than 50 cases [7]. The time needed to construct 
a boundary element model is about 10 — 15 minutes. 


































Figure 3: A software for registration. 


4 Discussion 

A software package was realized to reconstruct 
boundary element models from MR images for bio- 
electromagnetic inverse problems. The software has 
proved its usefulness in practice in many studies. 

A major problem with the MR imaging is that the 
ECG and MCG recordings and anatomic images are 
not acquired during the same session. Hence the pose 
of the patient may differ in separate measurements. 
The RMS error of registration is an indication of the 
problem. The error is not high if various error sources 
in the thorax region are considered, such as breathing, 
reattachement of the markers between MCG and MR 
measurements, shape changes between the measure¬ 
ments because of flexibility of the shoulders and the 
skin, and because of different shapes of the measure¬ 
ment beds (the MCG bed is flat while MRI bed is con¬ 
cave), and localization accuracy of the markers from 
MR images. However, in the scope of the MCG lo¬ 
calization accuracy, the registration error is relatively 
high and should be reduced; usually about 1 cm lo¬ 
calization accuracy of electric activity is required for 
clinical applicability in cardiac studies. A study for 
further reduction of different error sources is still on¬ 
going. The problems related to the elasticity of the 
head are smaller, because the skull of the human head 
is a rigid object. 
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